SUBROUTINE FFT_INTELL(nx,ny,init_free,dx,array2d,xfft2)

use MKL_DFTI
use VARKINDS
use FFT_INTEL, only : &
& ipar, dpar, Desc_Handle, xfft, &
& FFT_INTEL_ERROR_HANDLER

implicit none

! Input variables
integer(4), intent(in) :: nx, ny
integer(4), intent(inout) :: init_free

real(kind=kindreal), intent(in) :: dx
real(kind=kindreal), intent(inout) :: array2d(1:nx,1:ny)
real(kind=kindreal), intent(out) :: xfft2(1:nx)

! Local variables
integer(4) :: istatus
integer(4) :: i, j, k ! Loop indices

real(kind=kindreal), save :: pi
real(kind=kindreal) :: dk, work

!type(DFTI_DESCRIPTOR) :: desc_handle

!integer(4) :: precision, forward_domain
!integer(4) :: num_threads ! Number of OpenMP threads per descriptor
!integer(4), save :: fft_dimension = 1
!integer(4) :: i
!integer(4) :: jx, jy ! Sizes in x and y directions of a 3D array being
                     ! transformed in z direction

!real(8), allocatable :: x(:), x0(:)
!real(8) :: dscale

!allocate (x(1:nx), x0(1:nx))


!jx = 1
!jy = 1
!Status = DftiSetValue( Desc_Handle, DFTI_NUMBER_OF_TRANSFORMS, jx*jy )
!num_threads = 1
!Status = DftiSetValue( Desc_Handle, DFTI_NUMBER_OF_USER_THREADS, num_threads )
!Status = DftiSetValue( Desc_Handle, DFTI_INPUT_DISTANCE, nx)
!Status = DftiSetValue( Desc_Handle, DFTI_OUTPUT_DISTANCE, nx)
!Status = DftiSetValue( Desc_Handle, DFTI_FORWARD_SCALE, 1.d0 )
!Status = DftiSetValue( Desc_Handle, DFTI_BACKWARD_SCALE, 1.d0/float(nx) )
!Status = DftiGetValue( Desc_Handle, DFTI_FORWARD_SCALE, dscale )
!write (*,*) dscale
!Status = DftiGetValue( Desc_Handle, DFTI_BACKWARD_SCALE, dscale )
!write (*,*) dscale

!do i = 1, nx-2
!  x(i) = float(i)/float(nx)
!  x0(i) = x(i)
!enddo
!x(nx-1) = 0.d0
!x0(nx-1) = x(nx-1)

if (init_free == 0) then
! Initializing staggered cosine transform
!  allocate (dpar(1:5*(nx-1)/2+2))
  pi = 4.*datan(1.d0)
  allocate (xfft(1:nx))
  if (kindreal == 4) then
    i = DFTI_SINGLE
  elseif (kindreal == 8) then
    i = DFTI_DOUBLE
  endif
!  istatus = DftiCreateDescriptor( Desc_Handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, nx )
  call FFT_INTEL_ERROR_HANDLER('FFT_INTELL:DftiCreateDescriptor', &
  & DftiCreateDescriptor( Desc_Handle, i, DFTI_COMPLEX, 1, nx ) )
!  istatus = DftiCommitDescriptor( Desc_Handle )
  call FFT_INTEL_ERROR_HANDLER('FFT_INTELL:DftiCommitDescriptor', &
  & DftiCommitDescriptor( Desc_Handle ) )
!  call D_INIT_TRIG_TRANSFORM (nx-1, MKL_STAGGERED_COSINE_TRANSFORM, ipar, dpar, status)
!  if (status /= 0) call FFT_INTEL_ERROR_HANDLER('D_INIT_TRIG_TRANSFORM',status)
!  call D_COMMIT_TRIG_TRANSFORM (xfft, Desc_Handle, ipar, dpar, status)
!  if (status /= 0) call FFT_INTEL_ERROR_HANDLER('D_COMMIT_TRIG_TRANSFORM',status)
  init_free = 1
endif

!do i = 1, nx
!  xfft(i) = float(i)/float(nx)
!enddo
if (init_free == 1) then
  xfft2(:) = 0.d0
  do j = 1, ny
    xfft(1:nx) = cmplx(array2d(1:nx,j),0.)
!    xfft2(1:nx) = xfft(1:nx)
    call FFT_INTEL_ERROR_HANDLER('FFT_INTELL:DftiComputeForward', &
    & DftiComputeForward( desc_handle, xfft ) )
!    call FFT_INTEL_ERROR_HANDLER('FFT_INTELL:DftiComputeBackward', &
!    & DftiComputeBackward( desc_handle, xfft ) )
!    write(*,*) xfft - xfft2
!    STOP
    
!    do i = 1, nx
!      work = 0.
!      do k = 0, nx-1
!        work = work + 1./dble(nx)* &
!        & (real(xfft(k+1))*dcos(2.d0*pi*(i-1)*k/float(nx)) - &
!        &  aimag(xfft(k+1))*dsin(2.d0*pi*(i-1)*k/float(nx)))
!      enddo
!      write(*,*) array2d(i,j) - work
!    enddo
!    write(*,*) sum(array2d(1:nx,j)*array2d(1:nx,j)) - &
!    & sum(real(xfft(1:nx))*real(xfft(1:nx)))/float(nx) - &
!    & sum(aimag(xfft(1:nx))*aimag(xfft(1:nx)))/float(nx)

    do i = 1, nx
      if (i == 1 .or. i == nx) then
        dk = pi/(float(nx)*dx)
      else
        dk = 2.*pi/(float(nx)*dx)
      endif
      xfft2(i) = xfft2(i) + &
      &          (real(xfft(i))*real(xfft(i)) + &
      &          aimag(xfft(i))*aimag(xfft(i))) / float(nx)/ &
      &          dk
    enddo
!    write(*,*) 'A', array2d(1:nx,j), 'B', xfft, 'C', xfft2
!    STOP
  enddo
  xfft2(:) = xfft2(:)/float(ny)
endif


!call D_FORWARD_TRIG_TRANSFORM (xfft, Desc_Handle, ipar, dpar, status)
!if (status /= 0) call FFT_INTEL_ERROR_HANDLER('D_FORWARD_TRIG_TRANSFORM',status)
!call D_BACKWARD_TRIG_TRANSFORM (xfft, Desc_Handle, ipar, dpar, status)
!if (status /= 0) call FFT_INTEL_ERROR_HANDLER('D_BACKWARD_TRIG_TRANSFORM',status)

if (init_free == 2) then
! Making the allocated space free        
!  call FREE_TRIG_TRANSFORM (Desc_Handle, ipar, status)
!  if (status /= 0) call FFT_INTEL_ERROR_HANDLER('FREE_TRIG_TRANSFORM',status)
  call FFT_INTEL_ERROR_HANDLER('FFT_INTELL:DftiFreeDescriptor', &
  & DftiFreeDescriptor( Desc_Handle ) )
!   istatus = DftiFreeDescriptor( Desc_Handle )
!  deallocate (dpar)
  deallocate (xfft)
endif

!write(*,*) x - x0

!deallocate (x, x0)

! Debug
!call STOP_MODEL_MPI


END SUBROUTINE FFT_INTELL
